library(zoo)
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(xts)
library(CASdatasets)
## Loading required package: survival
library(dplyr)
##
## ######################### Warning from 'xts' package ##########################
## # #
## # The dplyr lag() function breaks how base R's lag() function is supposed to #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
## # source() into this session won't work correctly. #
## # #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
## # dplyr from breaking base R's lag() function. #
## # #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
## # #
## ###############################################################################
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:xts':
##
## first, last
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
library(ggplot2)
library(scales)
library(forcats)
library(kableExtra)
##
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
##
## group_rows
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
library(survival)
library(effects)
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(MASS)
library(broom)
## Warning: package 'broom' was built under R version 4.5.2
library(ggplot2)
library(mgcv)
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
## This is mgcv 1.9-4. For overview type '?mgcv'.
library(stringr)
library(statmod)
library(tweedie)
library(openxlsx)
library(glmmTMB)
## Warning in check_dep_version(dep_pkg = "TMB"): package version mismatch:
## glmmTMB was built with TMB package version 1.9.17
## Current TMB package version is 1.9.19
## Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)
##
## Attaching package: 'glmmTMB'
## The following object is masked from 'package:statmod':
##
## tweedie
library(cplm)
## Loading required package: coda
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
## Loading required package: splines
##
## Attaching package: 'cplm'
## The following object is masked from 'package:glmmTMB':
##
## VarCorr
## The following object is masked from 'package:nlme':
##
## VarCorr
#install.packages('xts')
install.packages("CASdatasets", repos = "http://dutangc.free.fr/pub/RRepos/pub/", type="source")
Carga de los 2 datasets
# Tabla frecuencia
data(freMTPL2freq)
freq <- freMTPL2freq
# Tabla severidad
data(freMTPL2sev)
sev <- freMTPL2sev
Exploración inicial –>
summary(freq)
## IDpol ClaimNb Exposure VehPower
## 1 : 1 Min. : 0.000 Min. :0.002732 Min. : 4.000
## 3 : 1 1st Qu.: 0.000 1st Qu.:0.180000 1st Qu.: 5.000
## 5 : 1 Median : 0.000 Median :0.490000 Median : 6.000
## 10 : 1 Mean : 0.039 Mean :0.528743 Mean : 6.455
## 11 : 1 3rd Qu.: 0.000 3rd Qu.:0.990000 3rd Qu.: 7.000
## 13 : 1 Max. :16.000 Max. :2.010000 Max. :15.000
## (Other):677985
## VehAge DrivAge BonusMalus VehBrand
## Min. : 0.000 Min. : 18.0 Min. : 50.00 B12 :166022
## 1st Qu.: 2.000 1st Qu.: 34.0 1st Qu.: 50.00 B1 :162730
## Median : 6.000 Median : 44.0 Median : 50.00 B2 :159854
## Mean : 7.044 Mean : 45.5 Mean : 59.76 B3 : 53393
## 3rd Qu.: 11.000 3rd Qu.: 55.0 3rd Qu.: 64.00 B5 : 34752
## Max. :100.000 Max. :100.0 Max. :230.00 B6 : 28545
## (Other): 72695
## VehGas Area Density
## Diesel :332128 A:103952 Min. : 1
## Regular:345863 B: 75457 1st Qu.: 92
## C:191874 Median : 393
## D:151592 Mean : 1792
## E:137163 3rd Qu.: 1658
## F: 17953 Max. :27000
##
## Region
## Centre :160595
## Rhone-Alpes : 84748
## Provence-Alpes-Cotes-D'Azur: 79315
## Ile-de-France : 69789
## Bretagne : 42120
## Pays-de-la-Loire : 38747
## (Other) :202677
names(freq)
## [1] "IDpol" "ClaimNb" "Exposure" "VehPower" "VehAge"
## [6] "DrivAge" "BonusMalus" "VehBrand" "VehGas" "Area"
## [11] "Density" "Region"
dim(freq)
## [1] 677991 12
12 variaable y 677991 registros.
Renombramiento variables –>
freq <- freq %>%
rename(
id_poliza = IDpol,
exposicion = Exposure,
numero_siniestros = ClaimNb,
potencia = VehPower,
edad_conductor = DrivAge,
edad_vehiculo = VehAge,
bonus_malus = BonusMalus,
marca = VehBrand,
combustible = VehGas,
area = Area,
region = Region
)
Exploración inicial –>
summary(sev)
## IDpol ClaimAmount
## 2241683: 16 Min. : 1
## 3253234: 11 1st Qu.: 686
## 3254353: 11 Median : 1172
## 2248174: 9 Mean : 2266
## 2239279: 8 3rd Qu.: 1212
## 2216294: 6 Max. :4075401
## (Other):26383
names(sev)
## [1] "IDpol" "ClaimAmount"
dim(sev)
## [1] 26444 2
2 variables y 26444 registros
Renombramiento variables –>
sev <- sev %>%
rename(
id_poliza = IDpol,
coste_siniestro = ClaimAmount
)
colSums(is.na(freq))
## id_poliza numero_siniestros exposicion potencia
## 0 0 0 0
## edad_vehiculo edad_conductor bonus_malus marca
## 0 0 0 0
## combustible area Density region
## 0 0 0 0
colSums(is.na(sev))
## id_poliza coste_siniestro
## 0 0
Analisis valores duplicados –>
any(duplicated(freq$IDpol))
## [1] FALSE
any(duplicated(sev$IDpol))
## [1] FALSE
sev <- sev %>%
group_by(id_poliza) %>%
summarise(
total_importe_siniestros = sum(coste_siniestro, na.rm = TRUE)
) %>%
ungroup()
Unificamos las dos tablas
polizas <- freq %>%
left_join(sev, by = "id_poliza")
Remplazamos los NA en el importe total de siniestros y en el número de siniestros por cero
polizas <- polizas %>%
mutate(
total_importe_siniestros = ifelse(is.na(total_importe_siniestros), 0, total_importe_siniestros),
)
Tasa de siniestralidad
Cuántos siniestros ocurren en promedio por póliza, ajustado por el tiempo que la póliza ha estado en vigor. Permite comprar pólizas que no han estado aseguradas el mismo tiempo.
\[ \text{tasa siniestros} = \frac{\text{numero siniestros}}{\text{exposición}} \]
polizas <- polizas %>%
mutate(tasa_siniestros = numero_siniestros / exposicion)
Prima pura estimada
Cantidad de dinero que, en promedio, se espera pagar en siniestros por póliza y unidad de exposición.
\[ \text{importe medio siniestro} = \begin{cases} \dfrac{\text{total importe siniestros}}{\text{numero siniestros}}, & \text{si } \text{numero siniestros} > 0 \\[6pt] 0, & \text{si } \text{numero siniestros} = 0 \end{cases} \]
\[ \text{prima pura estimada} = \text{tasa siniestros} \times \text{importe medio siniestro} \]
polizas <- polizas %>%
mutate(
importe_medio_siniestro = ifelse(numero_siniestros > 0, total_importe_siniestros / numero_siniestros, 0),
prima_pura_estimada = tasa_siniestros * importe_medio_siniestro
)
summary(polizas$exposicion)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.002732 0.180000 0.490000 0.528743 0.990000 2.010000
Distribución –>
ggplot(polizas, aes(x = exposicion)) +
geom_histogram(bins = 10, fill = "steelblue", color = "white") +
labs(title = "Distribución de la exposición", x = "Años asegurados", y = "Número de pólizas")+
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 12),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
ggplot(polizas, aes(x = "", y = exposicion)) +
geom_boxplot(fill = "steelblue", color = "grey2") +
labs(
title = "Distribución de la exposición",
y = "Exposición (años asegurados)",
x = ""
) +
theme_minimal()+
theme(
plot.title = element_text(hjust = 0.5, size = 12),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
Depuramos los valores superiores a un año –>
# Modificar las observaciones con exposición > 1, igualándolas a 1
polizas <- polizas %>%
mutate(exposicion = ifelse(exposicion > 1, 1, exposicion))
# Verificar el cambio
summary(polizas$exposicion)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.002732 0.180000 0.490000 0.528537 0.990000 1.000000
Distribución después de depurar la variable
ggplot(polizas, aes(x = exposicion)) +
geom_histogram(bins = 10, fill = "steelblue", color = "white") +
labs(title = "Distribución de la exposición", x = "Años asegurados", y = "Número de pólizas")+
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 12),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
Tabla descriptiva antes de la depuración
intervalos <- c(-Inf, 0, 1, 2, 3, 5, 10, 20, Inf)
etiquetas <- c("0", "1", "2", "3", "4-5", "6-10", "11-20", "Más de 20")
tabla_siniestros <- table(cut(polizas$numero_siniestros, breaks = intervalos, labels = etiquetas, right = TRUE)) %>%
as.data.frame()
colnames(tabla_siniestros) <- c("Intervalo", "Recuento")
tabla_siniestros$Porcentaje <- round(tabla_siniestros$Recuento / sum(tabla_siniestros$Recuento) * 100, 2)
tabla_siniestros
## Intervalo Recuento Porcentaje
## 1 0 653047 96.32
## 2 1 23571 3.48
## 3 2 1298 0.19
## 4 3 62 0.01
## 5 4-5 7 0.00
## 6 6-10 3 0.00
## 7 11-20 3 0.00
## 8 Más de 20 0 0.00
Eliminamos los valores superiores a 4 siniestros para evitar una sobre dispersión.
polizas <- polizas %>%
mutate(numero_siniestros = ifelse(numero_siniestros > 4, 4, numero_siniestros))
summary(polizas$numero_siniestros)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.03895 0.00000 4.00000
ggplot(polizas, aes(x = factor(numero_siniestros))) +
geom_bar(fill = "steelblue", color = "white") +
labs(
title = "Número de siniestros por póliza",
x = "Número de siniestros",
y = "Frecuencia"
) +
theme_minimal(base_size = 10) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 12),
plot.subtitle = element_text(hjust = 0.5, size = 9)
)
tabla_sin <- polizas %>%
count(numero_siniestros) %>%
mutate(
prop = n / sum(n),
etiqueta = percent(prop, accuracy = 0.1, decimal.mark = ","),
numero_siniestros = factor(numero_siniestros)
)
ggplot(tabla_sin, aes(x = numero_siniestros, y = n)) +
geom_col(fill = "steelblue", color = "white") +
geom_text(aes(label = etiqueta), vjust = -0.3, size = 3) +
labs(
title = "Número de siniestros por póliza",
x = "Número de siniestros en la póliza",
y = "Número de pólizas"
) +
scale_y_continuous(
labels = label_number(big.mark = ".", decimal.mark = ","),
expand = expansion(mult = c(0, 0.05))
) +
theme_minimal(base_size = 10) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 12),
plot.subtitle = element_text(hjust = 0.5, size = 9)
)
#### Tasa de sineistros
Proporción de pólizas con tasa igual a cero
tabla_cero <- polizas %>%
mutate(tasa_pos = if_else(tasa_siniestros == 0, "tasa siniestros = 0", "tasa sineistros > 0 siniestros")) %>%
count(tasa_pos) %>%
mutate(porcentaje = 100 * n / sum(n))
ggplot(tabla_cero, aes(x = tasa_pos, y = n)) +
geom_col(fill = "steelblue", color = "white") +
geom_text(aes(label = paste0(round(porcentaje, 1), "%")),
vjust = -0.5, size = 3) +
labs(
title = "Pólizas con tasa de siniestros nula y positiva",
x = "Situación",
y = "Número de pólizas"
)+
theme_minimal(base_size = 10) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 12),
plot.subtitle = element_text(hjust = 0.5, size = 9)
)
Histograma solo para tasas positivas y recortado
tasa_pos <- polizas %>%
filter(tasa_siniestros > 0,
tasa_siniestros <= 5) # recorte a 5; ajusta si hace falta
ggplot(tasa_pos, aes(x = tasa_siniestros)) +
geom_histogram(bins = 50, fill = "steelblue", color = "white") +
labs(
title = "Distribución de la tasa de siniestros (pólizas con tasa > 0,\nrecorte en 5)",
x = "Tasa de siniestros",
y = "Número de pólizas"
) +
theme_minimal(base_size = 10) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 12),
plot.subtitle = element_text(hjust = 0.5, size = 9)
)
tasa_pos <- polizas %>%
filter(tasa_siniestros > 0,
tasa_siniestros <= quantile(tasa_siniestros[tasa_siniestros > 0], 0.99))
ggplot(tasa_pos, aes(x = tasa_siniestros, y = exposicion)) +
geom_point(alpha = 0.15, color = "steelblue", size = 1) +
labs(
title = "Dispersión de la tasa de siniestros frente a la exposición",
x = "Tasa de siniestros",
y = "Exposición"
) +
theme_minimal(base_size = 10) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 12),
plot.subtitle = element_text(hjust = 0.5, size = 9)
)
ggplot(polizas %>% filter(total_importe_siniestros > 0),
aes(x = total_importe_siniestros)) +
geom_histogram(bins = 50, fill = "steelblue", color = "white") +
scale_x_log10(
labels = comma_format(accuracy = 1)) +
labs(title = "Importe total de siniestros (log)", x = "€ (escala log)", y = "Frecuencia") +
theme_minimal(base_size = 10) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 12),
plot.subtitle = element_text(hjust = 0.5, size = 9)
)
####Importe medio de siniestros
ggplot(polizas %>% filter(importe_medio_siniestro > 0),
aes(x = importe_medio_siniestro)) +
geom_histogram(bins = 50, fill = "steelblue", color = "white") +
scale_x_log10(
labels = comma_format(accuracy = 1)) +
labs(title = "Importe medio de siniestro (log)", x = "€ (escala log)", y = "Frecuencia") +
theme_minimal(base_size = 10) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 12),
plot.subtitle = element_text(hjust = 0.5, size = 9)
)
####Prima pura estimada
options(scipen = 999) # evita usar notación científica
summary(polizas$prima_pura_estimada)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 0.0 0.0 383.3 0.0 18524548.0
polizas %>%
filter(prima_pura_estimada > 0) %>%
ggplot(aes(x = prima_pura_estimada)) +
geom_histogram(bins = 60, fill = "steelblue", color = "white") +
scale_x_log10(labels = label_number(big.mark = ".", decimal.mark = ",")) +
labs(title = "Prima pura estimada (>0, escala log10)", x = "€ (log10)", y = "Frecuencia") +
theme_minimal(base_size = 10) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 12),
plot.subtitle = element_text(hjust = 0.5, size = 9)
)
Distribución edad del conductor
ggplot(polizas, aes(x = edad_conductor)) +
geom_histogram(bins = 50, fill = "steelblue", color = "white") +
labs(title = "Distribución de la edad del conductor", x = "Edad", y = "Frecuencia")+
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 12),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
Estudiamos cómo categorizar los datos:
exposicion_por_edad <- polizas %>%
group_by(edad_conductor) %>%
summarise(exposicion_total = sum(exposicion, na.rm = TRUE))
ggplot(exposicion_por_edad, aes(x = edad_conductor, y = exposicion_total)) +
geom_col(fill = "steelblue", color = "white") +
labs(
title = "Distribución de la exposición según la edad del conductor",
x = "Edad del conductor",
y = "Exposición"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 12, face = "bold"),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
siniestros_por_edad <- polizas %>%
group_by(edad_conductor) %>%
summarise(media_siniestros = mean(numero_siniestros, na.rm = TRUE))
ggplot(siniestros_por_edad, aes(x = edad_conductor, y = media_siniestros)) +
geom_col(fill = "steelblue", color = "white") +
labs(
title = "Número medio de siniestros por edad del conductor",
x = "Edad del conductor",
y = "Número medio de siniestros"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 12, face = "bold"),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
Categorizamos la variable:
polizas <- polizas %>%
mutate(
edad_conductor_cat = cut(
edad_conductor,
breaks = c(-Inf, 25, 40, 60, Inf),
labels = c("Jóvenes", "Adultos jóvenes", "Adultos","Mayores"),
right = TRUE
)
)
Estudiamos la variable edad categorizada:
porcentajes_edad <- polizas %>%
count(edad_conductor_cat) %>%
mutate(porcentaje = round(100 * n / sum(n), 1))
ggplot(porcentajes_edad, aes(x = edad_conductor_cat, y = n)) +
geom_col(fill = "steelblue", color = "white") +
geom_text(aes(label = paste0(porcentaje, "%")),
vjust = -0.5, size = 3) +
labs(
title = "Distribución por categoría de edad del conductor",
x = "Edad del conductor",
y = "Número de pólizas"
) +
expand_limits(y = max(porcentajes_edad$n) * 1.1) + # aumenta el eje Y un 10%
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(hjust = 0.5),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
Tasa media de siniestros por tramo de edad:
importe_por_edad <- polizas %>%
group_by(edad_conductor_cat) %>%
summarise(importe_medio = mean(importe_medio_siniestro, na.rm = TRUE)) %>%
ungroup()
ggplot(importe_por_edad,
aes(x = edad_conductor_cat, y = importe_medio)) +
geom_col(fill = "steelblue", color = "white") +
geom_text(aes(label = sprintf("%.0f €", importe_medio)),
vjust = -0.5, size = 3) +
labs(
title = "Importe medio del siniestro por grupo de edad del conductor",
x = "Grupo de edad del conductor",
y = "Importe medio (€)"
) +
expand_limits(y = max(importe_por_edad$importe_medio) * 1.15) + # ← aumenta eje Y un 15%
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(hjust = 0.5),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8),
panel.grid.minor = element_blank()
)
Depuramos la base de datos:
# Modificar las observaciones con exposición > 1, igualándolas a 1
polizas <- polizas %>%
mutate(edad_vehiculo = ifelse(edad_vehiculo > 30, 30, edad_vehiculo))
# Verificar el cambio
summary(polizas$edad_vehiculo)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 2.000 6.000 7.029 11.000 30.000
Distribución de la variable edad del vehículo depurada:
ggplot(polizas, aes(x = edad_vehiculo)) +
geom_histogram(bins = 30, fill = "steelblue", color = "white") +
labs(
title = "Distribución de la edad del vehículo",
x = "Años del vehículo",
y = "Número de pólizas"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 12),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
Categorizamos la variable
polizas <- polizas %>%
mutate(
edad_vehiculo_cat = case_when(
edad_vehiculo <= 5 ~ "Nuevo",
edad_vehiculo <= 12 ~ "Medio",
TRUE ~ "Viejo"
)
)
polizas$edad_vehiculo_cat <- factor(
polizas$edad_vehiculo_cat,
levels = c("Nuevo", "Medio", "Viejo")
)
Distribución variable edad del vehículo categorizada:
polizas %>%
count(edad_vehiculo_cat) %>%
mutate(porcentaje = 100 * n / sum(n)) %>%
ggplot(aes(x = edad_vehiculo_cat, y = porcentaje)) +
geom_col(fill = "steelblue") +
geom_text(aes(label = sprintf("%.1f%%", porcentaje)),
vjust = -0.3, size = 3) +
labs(
title = "Distribución porcentual por categorías de edad del vehículo",
x = "Categoría de edad del vehículo",
y = "Porcentaje (%)"
) +
theme_minimal(base_size = 10)
Tasa media de siniestros por categoría de edad del vehículo:
tasa_cat <- polizas %>%
group_by(edad_vehiculo_cat) %>%
summarise(
tasa_media = mean(tasa_siniestros, na.rm = TRUE),
.groups = "drop"
)
ggplot(tasa_cat, aes(x = edad_vehiculo_cat, y = tasa_media)) +
geom_col(fill = "steelblue") +
geom_text(aes(label = round(tasa_media, 3)),
vjust = -0.3, size = 3) +
labs(
title = "Tasa media de siniestros según la edad del vehículo",
x = "Categoría de edad del vehículo",
y = "Tasa media de siniestros"
)+
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 12),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
Importe medio de siniestro por categoría del vehículo
importe_cat <- polizas %>%
group_by(edad_vehiculo_cat) %>%
summarise(
importe_medio = mean(importe_medio_siniestro, na.rm = TRUE),
.groups = "drop"
)
ggplot(importe_cat, aes(x = edad_vehiculo_cat, y = importe_medio)) +
geom_col(fill = "steelblue") +
geom_text(aes(label = paste0(round(importe_medio, 0), " €")),
vjust = -0.3, size = 3) +
labs(
title = "Importe medio de siniestro según la edad del vehículo",
x = "Categoría de edad del vehículo",
y = "Importe medio (€)"
)+
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 12),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8))
ggplot(polizas, aes(x = potencia)) +
geom_bar(fill = "steelblue", color = "white") +
labs(
title = "Distribución de la potencia del vehículo",
x = "Potencia",
y = "Número de pólizas"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 12),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
siniestros_medios_por_pot <- polizas %>%
group_by(potencia) %>%
summarise(
siniestros_medios = mean(numero_siniestros, na.rm = TRUE)
) %>%
ungroup()
ggplot(siniestros_medios_por_pot, aes(x = potencia, y = siniestros_medios)) +
geom_col(fill = "steelblue", color = "white") +
geom_text(aes(label = sprintf("%.3f", siniestros_medios)),
vjust = -0.5, size = 3) +
labs(
title = "Número medio de siniestros por potencia del vehículo",
x = "Potencia del vehículo",
y = "Número medio de siniestros por póliza"
)+
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 12),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
tasa_por_pot <- polizas %>%
group_by(potencia) %>%
summarise(
tasa_media = sum(numero_siniestros, na.rm = TRUE) / sum(exposicion, na.rm = TRUE)
) %>%
ungroup()
ggplot(tasa_por_pot, aes(x = potencia, y = tasa_media)) +
geom_line(color = "steelblue", linewidth = 1.2) +
geom_point(color = "steelblue", size = 2) +
geom_text(aes(label = sprintf("%.3f", tasa_media)),
vjust = -0.8, size = 3) +
labs(
title = "Tasa media de siniestralidad por potencia del vehículo",
x = "Potencia del vehículo",
y = "Tasa media de siniestralidad"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 12),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
Potencia en grupo
polizas <- polizas %>%
mutate(
potencia_grupo = case_when(
potencia <= 5 ~ "Baja",
potencia > 5 & potencia <= 9 ~ "Media",
potencia > 9 & potencia <= 12 ~ "Alta",
potencia > 12 ~ "Muy Alta",
TRUE ~ NA_character_
)
)
ggplot(polizas, aes(x = potencia_grupo)) +
geom_bar(fill = "steelblue", color = "white") +
geom_text(stat = "count", aes(label = scales::comma(..count..)),
vjust = -0.4, size = 2.8) + # añade etiquetas encima de las barras
labs(
title = "Distribución por grupo de potencia",
x = "Grupo de potencia",
y = "Número de pólizas"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 12),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
### Modelización de la frecuencia
Preparamos y comprobamos la variable respuesta
summary(polizas$numero_siniestros)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.03895 0.00000 4.00000
table(polizas$numero_siniestros)
##
## 0 1 2 3 4
## 653047 23571 1298 62 13
mean(polizas$numero_siniestros)
## [1] 0.03894594
var(polizas$numero_siniestros)
## [1] 0.04203695
mean(polizas$numero_siniestros == 0)
## [1] 0.963209
Probamos la distribución Poisson
Función auxiliar para la sobredispersión:
dispersion <- function(model) {
sum(residuals(model, type = "pearson")^2) / df.residual(model)
}
polizas$potencia_grupo <- factor(polizas$potencia_grupo, ordered = FALSE)
Ajustamos el modelo Poisson básico
modelo_posion_base <- glm(numero_siniestros ~ edad_conductor_cat + edad_vehiculo_cat + potencia_grupo + combustible + bonus_malus,
offset = log(exposicion),
family = poisson,
data = polizas)
summary(modelo_posion_base)
##
## Call:
## glm(formula = numero_siniestros ~ edad_conductor_cat + edad_vehiculo_cat +
## potencia_grupo + combustible + bonus_malus, family = poisson,
## data = polizas, offset = log(exposicion))
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) -4.1510498 0.0441119 -94.103
## edad_conductor_catAdultos jóvenes -0.0833998 0.0245830 -3.393
## edad_conductor_catAdultos 0.2175425 0.0256199 8.491
## edad_conductor_catMayores 0.1178927 0.0290862 4.053
## edad_vehiculo_catMedio 0.0521224 0.0136185 3.827
## edad_vehiculo_catViejo -0.2121828 0.0177988 -11.921
## potencia_grupoBaja -0.1932816 0.0238202 -8.114
## potencia_grupoMedia -0.1261790 0.0227993 -5.534
## potencia_grupoMuy Alta 0.0859290 0.0586893 1.464
## combustibleRegular -0.1313401 0.0126962 -10.345
## bonus_malus 0.0274305 0.0003405 80.566
## Pr(>|z|)
## (Intercept) < 0.0000000000000002 ***
## edad_conductor_catAdultos jóvenes 0.000692 ***
## edad_conductor_catAdultos < 0.0000000000000002 ***
## edad_conductor_catMayores 0.000050518233659135 ***
## edad_vehiculo_catMedio 0.000130 ***
## edad_vehiculo_catViejo < 0.0000000000000002 ***
## potencia_grupoBaja 0.000000000000000489 ***
## potencia_grupoMedia 0.000000031238680864 ***
## potencia_grupoMuy Alta 0.143158
## combustibleRegular < 0.0000000000000002 ***
## bonus_malus < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 171311 on 677990 degrees of freedom
## Residual deviance: 164937 on 677980 degrees of freedom
## AIC: 215722
##
## Number of Fisher Scoring iterations: 6
dispersion(modelo_posion_base)
## [1] 1.664483
AIC(modelo_posion_base)
## [1] 215721.7
par(mfrow = c(2, 2)); plot(modelo_posion_base); par(mfrow = c(1, 1))
Probamos la distribución Binomial negativa
modelo_binomial_negativo <- glm.nb(
numero_siniestros ~ edad_conductor_cat +
edad_vehiculo_cat +
potencia_grupo +
combustible +
bonus_malus +
offset(log(exposicion)),
data = polizas
)
summary(modelo_binomial_negativo)
##
## Call:
## glm.nb(formula = numero_siniestros ~ edad_conductor_cat + edad_vehiculo_cat +
## potencia_grupo + combustible + bonus_malus + offset(log(exposicion)),
## data = polizas, init.theta = 1.043013931, link = log)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) -4.1850273 0.0468518 -89.325
## edad_conductor_catAdultos jóvenes -0.0778844 0.0257736 -3.022
## edad_conductor_catAdultos 0.2290706 0.0270195 8.478
## edad_conductor_catMayores 0.1264174 0.0305558 4.137
## edad_vehiculo_catMedio 0.0539929 0.0140405 3.846
## edad_vehiculo_catViejo -0.2128577 0.0182840 -11.642
## potencia_grupoBaja -0.1907287 0.0245545 -7.768
## potencia_grupoMedia -0.1248341 0.0235090 -5.310
## potencia_grupoMuy Alta 0.0892657 0.0604698 1.476
## combustibleRegular -0.1320668 0.0130835 -10.094
## bonus_malus 0.0278758 0.0003724 74.851
## Pr(>|z|)
## (Intercept) < 0.0000000000000002 ***
## edad_conductor_catAdultos jóvenes 0.00251 **
## edad_conductor_catAdultos < 0.0000000000000002 ***
## edad_conductor_catMayores 0.000035146351850 ***
## edad_vehiculo_catMedio 0.00012 ***
## edad_vehiculo_catViejo < 0.0000000000000002 ***
## potencia_grupoBaja 0.000000000000008 ***
## potencia_grupoMedia 0.000000109589658 ***
## potencia_grupoMuy Alta 0.13989
## combustibleRegular < 0.0000000000000002 ***
## bonus_malus < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.043) family taken to be 1)
##
## Null deviance: 151138 on 677990 degrees of freedom
## Residual deviance: 145155 on 677980 degrees of freedom
## AIC: 215185
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.0430
## Std. Err.: 0.0576
##
## 2 x log-likelihood: -215161.2840
AIC(modelo_binomial_negativo)
## [1] 215185.3
dispersion(modelo_binomial_negativo)
## [1] 1.622827
par(mfrow = c(2, 2)); plot(modelo_binomial_negativo); par(mfrow = c(1, 1))
Estructura del predictor usando Poisson
Edad del conductor
m_pois_age_lin <- glm(
numero_siniestros ~ edad_conductor +
edad_vehiculo_cat +
potencia_grupo +
combustible +
bonus_malus +
offset(log(exposicion)),
family = poisson,
data = polizas
)
AIC(modelo_posion_base, m_pois_age_lin)
## df AIC
## modelo_posion_base 11 215721.7
## m_pois_age_lin 9 215985.6
dispersion(m_pois_age_lin)
## [1] 1.67128
m_pois_age_quad <- glm(
numero_siniestros ~ edad_conductor + I(edad_conductor^2) +
edad_vehiculo_cat +
potencia_grupo +
combustible +
bonus_malus +
offset(log(exposicion)),
family = poisson,
data = polizas
)
AIC(modelo_posion_base, m_pois_age_quad)
## df AIC
## modelo_posion_base 11 215721.7
## m_pois_age_quad 10 215890.0
dispersion(m_pois_age_quad)
## [1] 1.665208
library(splines)
m_pois_age_spline <- glm(
numero_siniestros ~ ns(edad_conductor, df = 4) +
edad_vehiculo_cat +
potencia_grupo +
combustible +
bonus_malus +
offset(log(exposicion)),
family = poisson,
data = polizas
)
AIC(modelo_posion_base, m_pois_age_spline)
## df AIC
## modelo_posion_base 11 215721.7
## m_pois_age_spline 12 215544.5
dispersion(m_pois_age_spline)
## [1] 1.661709
Edad vehículo
m_pois_veh_lin <- glm(
numero_siniestros ~ edad_conductor_cat + # o tu forma elegida arriba
edad_vehiculo +
potencia_grupo +
combustible +
bonus_malus +
offset(log(exposicion)),
family = poisson,
data = polizas
)
AIC(m_pois_veh_lin)
## [1] 215840.3
dispersion(m_pois_veh_lin)
## [1] 1.659176
m_pois_veh_quad <- glm(
numero_siniestros ~ edad_conductor_cat +
edad_vehiculo + I(edad_vehiculo^2) +
potencia_grupo +
combustible +
bonus_malus +
offset(log(exposicion)),
family = poisson,
data = polizas
)
AIC(m_pois_veh_quad)
## [1] 215676.7
dispersion(m_pois_veh_quad)
## [1] 1.662875
m_pois_veh_spline <- glm(
numero_siniestros ~ edad_conductor_cat +
ns(edad_vehiculo, df = 4) +
potencia_grupo +
combustible +
bonus_malus +
offset(log(exposicion)),
family = poisson,
data = polizas
)
AIC(m_pois_veh_spline)
## [1] 215661
dispersion(m_pois_veh_spline)
## [1] 1.66449
Potencia
m_pois_pot_lin <- glm(
numero_siniestros ~ edad_conductor_cat +
edad_vehiculo_cat +
potencia +
combustible +
bonus_malus +
offset(log(exposicion)),
family = poisson,
data = polizas
)
AIC(m_pois_pot_lin)
## [1] 215696.6
dispersion(m_pois_pot_lin)
## [1] 1.665391
m_pois_veh_quad_lin <- glm(
numero_siniestros ~ edad_conductor_cat +
edad_vehiculo + I(edad_vehiculo^2) +
potencia +
combustible +
bonus_malus +
offset(log(exposicion)),
family = poisson,
data = polizas
)
AIC(m_pois_veh_quad_lin)
## [1] 215648.9
dispersion(m_pois_veh_quad_lin)
## [1] 1.664081
m_pois_pot_quad <- glm(
numero_siniestros ~ edad_conductor_cat +
edad_vehiculo_cat +
potencia + I(potencia^2) +
combustible +
bonus_malus +
offset(log(exposicion)),
family = poisson,
data = polizas
)
AIC(m_pois_pot_quad)
## [1] 215698.1
dispersion(m_pois_pot_quad)
## [1] 1.665872
m_pois_pot_log <- glm(
numero_siniestros ~ edad_conductor_cat +
edad_vehiculo_cat +
log(potencia) +
combustible +
bonus_malus +
offset(log(exposicion)),
family = poisson,
data = polizas
)
AIC(m_pois_pot_log)
## [1] 215696.9
dispersion(m_pois_pot_log)
## [1] 1.666867
Actualizo modelo base
modelo_poison_base_1 <- glm(
numero_siniestros ~ edad_conductor_cat +
edad_vehiculo + I(edad_vehiculo^2) +
potencia +
combustible +
bonus_malus +
offset(log(exposicion)),
family = poisson,
data = polizas
)
summary(modelo_poison_base_1)
##
## Call:
## glm(formula = numero_siniestros ~ edad_conductor_cat + edad_vehiculo +
## I(edad_vehiculo^2) + potencia + combustible + bonus_malus +
## offset(log(exposicion)), family = poisson, data = polizas)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) -4.5529118 0.0454218 -100.236
## edad_conductor_catAdultos jóvenes -0.0852246 0.0245875 -3.466
## edad_conductor_catAdultos 0.2150798 0.0256282 8.392
## edad_conductor_catMayores 0.1201094 0.0290828 4.130
## edad_vehiculo 0.0323555 0.0036902 8.768
## I(edad_vehiculo^2) -0.0025241 0.0002015 -12.529
## potencia 0.0338055 0.0030717 11.005
## combustibleRegular -0.1336424 0.0124918 -10.698
## bonus_malus 0.0274244 0.0003404 80.567
## Pr(>|z|)
## (Intercept) < 0.0000000000000002 ***
## edad_conductor_catAdultos jóvenes 0.000528 ***
## edad_conductor_catAdultos < 0.0000000000000002 ***
## edad_conductor_catMayores 0.0000363 ***
## edad_vehiculo < 0.0000000000000002 ***
## I(edad_vehiculo^2) < 0.0000000000000002 ***
## potencia < 0.0000000000000002 ***
## combustibleRegular < 0.0000000000000002 ***
## bonus_malus < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 171311 on 677990 degrees of freedom
## Residual deviance: 164868 on 677982 degrees of freedom
## AIC: 215649
##
## Number of Fisher Scoring iterations: 6
AIC(modelo_poison_base_1)
## [1] 215648.9
dispersion(modelo_poison_base_1)
## [1] 1.664081
Combustible
summary(modelo_poison_base_1)$coefficients["combustibleRegular",]
## Estimate Std. Error
## -0.13364239045306164355153555334255 0.01249178909848610782851174860753
## z value Pr(>|z|)
## -10.69841872924815007195320504251868 0.00000000000000000000000001035298
m_pois_no_comb <- glm(
numero_siniestros ~ edad_conductor_cat +
edad_vehiculo_cat +
potencia +
bonus_malus +
offset(log(exposicion)),
family = poisson,
data = polizas
)
AIC(modelo_poison_base_1, m_pois_no_comb)
## df AIC
## modelo_poison_base_1 9 215648.9
## m_pois_no_comb 8 215808.6
dispersion(m_pois_no_comb)
## [1] 1.65824
anova(modelo_poison_base_1, m_pois_no_comb, test="Chisq")
## Analysis of Deviance Table
##
## Model 1: numero_siniestros ~ edad_conductor_cat + edad_vehiculo + I(edad_vehiculo^2) +
## potencia + combustible + bonus_malus + offset(log(exposicion))
## Model 2: numero_siniestros ~ edad_conductor_cat + edad_vehiculo_cat +
## potencia + bonus_malus + offset(log(exposicion))
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 677982 164868
## 2 677983 165030 -1 -161.71 < 0.00000000000000022 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m_pois_interact <- glm(
numero_siniestros ~ edad_conductor_cat +
edad_vehiculo + I(edad_vehiculo^2)+
potencia*combustible +
bonus_malus +
offset(log(exposicion)),
family = poisson,
data = polizas
)
AIC(modelo_poison_base_1, m_pois_interact)
## df AIC
## modelo_poison_base_1 9 215648.9
## m_pois_interact 10 215635.3
dispersion(m_pois_interact)
## [1] 1.66635
anova(modelo_poison_base_1, m_pois_interact, test="Chisq")
## Analysis of Deviance Table
##
## Model 1: numero_siniestros ~ edad_conductor_cat + edad_vehiculo + I(edad_vehiculo^2) +
## potencia + combustible + bonus_malus + offset(log(exposicion))
## Model 2: numero_siniestros ~ edad_conductor_cat + edad_vehiculo + I(edad_vehiculo^2) +
## potencia * combustible + bonus_malus + offset(log(exposicion))
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 677982 164868
## 2 677981 164853 1 15.562 0.00007983 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Actualizo modelo base
m_pois_final <- glm(
numero_siniestros ~ edad_conductor_cat +
edad_vehiculo + I(edad_vehiculo^2)+
potencia*combustible +
bonus_malus +
offset(log(exposicion)),
family = poisson,
data = polizas
)
summary(m_pois_final)
##
## Call:
## glm(formula = numero_siniestros ~ edad_conductor_cat + edad_vehiculo +
## I(edad_vehiculo^2) + potencia * combustible + bonus_malus +
## offset(log(exposicion)), family = poisson, data = polizas)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) -4.4547759 0.0518204 -85.966
## edad_conductor_catAdultos jóvenes -0.0839910 0.0245884 -3.416
## edad_conductor_catAdultos 0.2192898 0.0256506 8.549
## edad_conductor_catMayores 0.1237228 0.0290992 4.252
## edad_vehiculo 0.0313359 0.0036940 8.483
## I(edad_vehiculo^2) -0.0024707 0.0002016 -12.258
## potencia 0.0189153 0.0049032 3.858
## combustibleRegular -0.2935956 0.0424651 -6.914
## bonus_malus 0.0274189 0.0003404 80.543
## potencia:combustibleRegular 0.0246730 0.0062642 3.939
## Pr(>|z|)
## (Intercept) < 0.0000000000000002 ***
## edad_conductor_catAdultos jóvenes 0.000636 ***
## edad_conductor_catAdultos < 0.0000000000000002 ***
## edad_conductor_catMayores 0.00002120980838 ***
## edad_vehiculo < 0.0000000000000002 ***
## I(edad_vehiculo^2) < 0.0000000000000002 ***
## potencia 0.000114 ***
## combustibleRegular 0.00000000000472 ***
## bonus_malus < 0.0000000000000002 ***
## potencia:combustibleRegular 0.00008190507532 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 171311 on 677990 degrees of freedom
## Residual deviance: 164853 on 677981 degrees of freedom
## AIC: 215635
##
## Number of Fisher Scoring iterations: 6
AIC(m_pois_final)
## [1] 215635.3
dispersion(m_pois_final)
## [1] 1.66635
par(mfrow=c(2,2)); plot(m_pois_final); par(mfrow=c(1,1))
Ajustamos Binomial negativa
m_nb_final <- glm.nb(
numero_siniestros ~ edad_conductor_cat +
edad_vehiculo + I(edad_vehiculo^2)+
potencia*combustible +
bonus_malus +
offset(log(exposicion)),
data = polizas
)
summary(m_nb_final)
##
## Call:
## glm.nb(formula = numero_siniestros ~ edad_conductor_cat + edad_vehiculo +
## I(edad_vehiculo^2) + potencia * combustible + bonus_malus +
## offset(log(exposicion)), data = polizas, init.theta = 1.048476778,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) -4.4838603 0.0546565 -82.037
## edad_conductor_catAdultos jóvenes -0.0783767 0.0257741 -3.041
## edad_conductor_catAdultos 0.2307577 0.0270422 8.533
## edad_conductor_catMayores 0.1325212 0.0305647 4.336
## edad_vehiculo 0.0315287 0.0037834 8.334
## I(edad_vehiculo^2) -0.0024803 0.0002061 -12.036
## potencia 0.0184044 0.0050530 3.642
## combustibleRegular -0.2969096 0.0437502 -6.786
## bonus_malus 0.0278638 0.0003723 74.845
## potencia:combustibleRegular 0.0250866 0.0064541 3.887
## Pr(>|z|)
## (Intercept) < 0.0000000000000002 ***
## edad_conductor_catAdultos jóvenes 0.002359 **
## edad_conductor_catAdultos < 0.0000000000000002 ***
## edad_conductor_catMayores 0.0000145261777 ***
## edad_vehiculo < 0.0000000000000002 ***
## I(edad_vehiculo^2) < 0.0000000000000002 ***
## potencia 0.000270 ***
## combustibleRegular 0.0000000000115 ***
## bonus_malus < 0.0000000000000002 ***
## potencia:combustibleRegular 0.000102 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.0485) family taken to be 1)
##
## Null deviance: 151219 on 677990 degrees of freedom
## Residual deviance: 145153 on 677981 degrees of freedom
## AIC: 215102
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.0485
## Std. Err.: 0.0581
##
## 2 x log-likelihood: -215080.3840
comparamos con Poisson
tabla1 <- AIC(m_pois_final, m_nb_final)
dispersion(m_nb_final)
## [1] 1.624849
Revisamos los residuos
par(mfrow=c(2,2)); plot(m_nb_final); par(mfrow=c(1,1))
Gráficos básicos de diagnóstico
par(mfrow = c(2, 2))
plot(m_nb_final)
par(mfrow = c(1, 1))
res_nb <- residuals(m_nb_final, type = "pearson")
lev_nb <- hatvalues(m_nb_final)
cook_nb <- cooks.distance(m_nb_final)
sort(cook_nb, decreasing=TRUE)[1:10]
## 259414 39811 141089 558610 574655 152687
## 0.0032881621 0.0031219588 0.0025989982 0.0014153098 0.0009295472 0.0009032286
## 469045 541664 589135 259341
## 0.0009012370 0.0008936455 0.0008918052 0.0008802308
eff_df <- as.data.frame(effect("edad_vehiculo", m_nb_final))
ggplot(eff_df, aes(x = edad_vehiculo, y = fit)) +
geom_line(fill = "steelblue",size = 1) +
geom_ribbon(fill = "steelblue",aes(ymin = lower, ymax = upper), alpha = 0.2) +
labs(
title = "Efecto no lineal de la edad del vehículo sobre la frecuencia siniestral",
x = "Edad del vehículo",
y = "Frecuencia esperada de siniestros"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 12),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning in geom_line(fill = "steelblue", size = 1): Ignoring unknown
## parameters: `fill`
Definir la variable severidad y conjunto de datos
datos_sev <- polizas %>%
filter(numero_siniestros > 0) %>%
mutate(
sev_media = total_importe_siniestros / numero_siniestros
)
summary(datos_sev$sev_media)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.0 710.6 1172.0 2222.3 1228.1 4075400.6
quantile(datos_sev$sev_media, probs = c(0.90, 0.95, 0.99, 0.995, 0.999), na.rm = TRUE)
## 90% 95% 99% 99.5% 99.9%
## 2742.509 4659.463 16327.001 34564.646 133304.840
corte <- quantile(datos_sev$sev_media, probs = 0.995, na.rm = TRUE)
datos_sev$sev_trunc <- pmin(datos_sev$sev_media, corte)
summary(datos_sev$sev_trunc)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.0 710.6 1172.0 1683.8 1228.1 34564.6
prop_trunc <- mean(datos_sev$sev_media > corte)
prop_trunc
## [1] 0.005011225
Exploratorio rápido de la severidad (para justificar Gamma/log)
ggplot(datos_sev, aes(x = sev_trunc)) +
geom_histogram(bins = 50, fill = "#4C72B0", alpha = 0.9) +
labs(
title = "Severidad truncada",
x = "Severidad",
y = "Número de pólizas"
) + theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 12),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
ggplot(datos_sev, aes(x = log(sev_trunc))) +
geom_histogram(bins = 50, fill = "#4C72B0", alpha = 0.9) +
labs(
title = "Log-severidad",
x = "log(Severidad)",
y = "Número de pólizas"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 12),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
Definimos los factores explicativos
datos_sev <- datos_sev %>%
mutate(
edad_conductor_cat = factor(edad_conductor_cat),
edad_vehiculo_cat = factor(edad_vehiculo_cat),
potencia_grupo = factor(potencia_grupo),
combustible = factor(combustible),
bonus_malus_num = as.numeric(as.character(bonus_malus))
)
Selección del predictor lineal
mod_sev_0 <- glm(
sev_trunc ~ edad_conductor_cat +
edad_vehiculo_cat +
potencia_grupo +
bonus_malus,
family = Gamma(link = "log"),
data = datos_sev,
weights = numero_siniestros
)
summary(mod_sev_0)
##
## Call:
## glm(formula = sev_trunc ~ edad_conductor_cat + edad_vehiculo_cat +
## potencia_grupo + bonus_malus, family = Gamma(link = "log"),
## data = datos_sev, weights = numero_siniestros)
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 7.4775042 0.0900126 83.072
## edad_conductor_catAdultos jóvenes -0.0846146 0.0490522 -1.725
## edad_conductor_catAdultos -0.0732701 0.0516660 -1.418
## edad_conductor_catMayores -0.0292565 0.0578605 -0.506
## edad_vehiculo_catMedio -0.0901533 0.0266370 -3.385
## edad_vehiculo_catViejo -0.1126957 0.0346615 -3.251
## potencia_grupoBaja -0.1150261 0.0462989 -2.484
## potencia_grupoMedia -0.0502134 0.0445056 -1.128
## potencia_grupoMuy Alta 0.0677571 0.1141975 0.593
## bonus_malus 0.0023140 0.0007252 3.191
## Pr(>|t|)
## (Intercept) < 0.0000000000000002 ***
## edad_conductor_catAdultos jóvenes 0.084541 .
## edad_conductor_catAdultos 0.156159
## edad_conductor_catMayores 0.613115
## edad_vehiculo_catMedio 0.000714 ***
## edad_vehiculo_catViejo 0.001150 **
## potencia_grupoBaja 0.012983 *
## potencia_grupoMedia 0.259225
## potencia_grupoMuy Alta 0.552964
## bonus_malus 0.001420 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 3.806349)
##
## Null deviance: 30402 on 24943 degrees of freedom
## Residual deviance: 30235 on 24934 degrees of freedom
## AIC: 446180
##
## Number of Fisher Scoring iterations: 6
Añadimos combustible
mod_sev_1 <- glm(
sev_trunc ~ edad_conductor_cat +
edad_vehiculo_cat +
potencia_grupo +
bonus_malus +
combustible,
family = Gamma(link = "log"),
data = datos_sev,
weights = numero_siniestros
)
summary(mod_sev_1)
##
## Call:
## glm(formula = sev_trunc ~ edad_conductor_cat + edad_vehiculo_cat +
## potencia_grupo + bonus_malus + combustible, family = Gamma(link = "log"),
## data = datos_sev, weights = numero_siniestros)
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 7.4696116 0.0904991 82.538
## edad_conductor_catAdultos jóvenes -0.0833713 0.0490685 -1.699
## edad_conductor_catAdultos -0.0727989 0.0516524 -1.409
## edad_conductor_catMayores -0.0309608 0.0578777 -0.535
## edad_vehiculo_catMedio -0.0910682 0.0266443 -3.418
## edad_vehiculo_catViejo -0.1164238 0.0349163 -3.334
## potencia_grupoBaja -0.1166356 0.0463453 -2.517
## potencia_grupoMedia -0.0489151 0.0445077 -1.099
## potencia_grupoMuy Alta 0.0599634 0.1144862 0.524
## bonus_malus 0.0022927 0.0007253 3.161
## combustibleRegular 0.0200524 0.0246449 0.814
## Pr(>|t|)
## (Intercept) < 0.0000000000000002 ***
## edad_conductor_catAdultos jóvenes 0.089316 .
## edad_conductor_catAdultos 0.158729
## edad_conductor_catMayores 0.592700
## edad_vehiculo_catMedio 0.000632 ***
## edad_vehiculo_catViejo 0.000856 ***
## potencia_grupoBaja 0.011853 *
## potencia_grupoMedia 0.271768
## potencia_grupoMuy Alta 0.600450
## bonus_malus 0.001575 **
## combustibleRegular 0.415853
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 3.803772)
##
## Null deviance: 30402 on 24943 degrees of freedom
## Residual deviance: 30232 on 24933 degrees of freedom
## AIC: 446179
##
## Number of Fisher Scoring iterations: 6
AIC(mod_sev_0, mod_sev_1)
## df AIC
## mod_sev_0 11 446179.9
## mod_sev_1 12 446179.3
Modelo final
mod_sev_final <- mod_sev_0
Diagnóstico del modelo final Revisamos la dispersión y calidad global del modelo
dispersion <- mod_sev_final$deviance / mod_sev_final$df.residual
dispersion
## [1] 1.212588
Interpretar los coeficientes
coef_sev <- coef(mod_sev_final)
relatividades <- exp(coef_sev)
round(relatividades, 3)
## (Intercept) edad_conductor_catAdultos jóvenes
## 1767.823 0.919
## edad_conductor_catAdultos edad_conductor_catMayores
## 0.929 0.971
## edad_vehiculo_catMedio edad_vehiculo_catViejo
## 0.914 0.893
## potencia_grupoBaja potencia_grupoMedia
## 0.891 0.951
## potencia_grupoMuy Alta bonus_malus
## 1.070 1.002
Diagnóstico del modelo de severidad
# Preparar datos de diagnóstico
datos_diag <- datos_sev %>%
mutate(
fitted = fitted(mod_sev_0),
resid_dev = residuals(mod_sev_0, type = "deviance"),
resid_pearson = residuals(mod_sev_0, type = "pearson"),
cook = cooks.distance(mod_sev_0)
)
# Residuos vs Ajustados
p1 <- ggplot(datos_diag, aes(x = fitted, y = resid_dev)) +
geom_point(alpha = 0.25) +
geom_hline(yintercept = 0, linetype = 2) +
geom_smooth(se = FALSE, method = "loess") +
labs(title = "Residuos devianza vs ajustados", x = "μ̂", y = "Residuos (dev)") +
theme_minimal()
# Q-Q plot residuos Pearson
p2 <- ggplot(datos_diag %>% filter(is.finite(resid_pearson)),
aes(sample = resid_pearson)) +
stat_qq(alpha = 0.25) +
stat_qq_line() +
labs(title = "Q-Q residuos Pearson", x = "Cuantiles teóricos", y = "Residuos") +
theme_minimal()
# Cook's distance
n <- nobs(mod_sev_0)
umbral <- 4 / n
p3 <- ggplot(datos_diag, aes(x = seq_along(cook), y = cook)) +
geom_col() +
geom_hline(yintercept = umbral, linetype = 2) +
labs(title = "Distancia de Cook", x = "Observación", y = "Cook") +
theme_minimal()
# Calibración por deciles (O/E)
datos_cal <- datos_sev %>%
mutate(
sev_esp = fitted(mod_sev_0),
sev_obs = sev_trunc,
w = numero_siniestros
) %>%
filter(is.finite(sev_esp), is.finite(sev_obs), w > 0) %>%
mutate(decil = ntile(sev_esp, 10)) %>%
group_by(decil) %>%
summarise(
sev_obs = sum(sev_obs * w) / sum(w),
sev_esp = sum(sev_esp * w) / sum(w),
OE = sev_obs / sev_esp,
.groups = "drop"
)
p4 <- ggplot(datos_cal, aes(x = decil, y = OE)) +
geom_line() +
geom_point(size = 2) +
geom_hline(yintercept = 1, linetype = 2) +
scale_x_continuous(breaks = 1:10) +
labs(title = "Calibración por deciles (O/E)", x = "Decil μ̂", y = "O/E") +
theme_minimal()
Asiganamos la variable respuesta
polizas <- polizas %>%
mutate(
S_bruto = ifelse(numero_siniestros == 0, 0,
total_importe_siniestros)
)
summary(polizas$S_bruto)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 0.0 0.0 88.4 0.0 4075400.6
mean(polizas$S_bruto == 0)
## [1] 0.963209
p_trunc <- 0.995
c_trunc <- quantile(
polizas$S_bruto[polizas$S_bruto > 0],
probs = p_trunc,
na.rm = TRUE
)
c_trunc
## 99.5%
## 35630.33
polizas <- polizas %>%
mutate(
S_trunc = pmin(S_bruto, as.numeric(c_trunc))
)
summary(polizas$S_trunc)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 0.00 0.00 66.13 0.00 35630.33
max(polizas$S_trunc)
## [1] 35630.33
c(
umbral_truncado = as.numeric(c_trunc),
proporcion_truncadas = mean(polizas$S_bruto > c_trunc),
impacto_medio = mean(polizas$S_bruto - polizas$S_trunc)
)
## umbral_truncado proporcion_truncadas impacto_medio
## 35630.3308000000 0.0001843682 22.2357667581
** Especificación y estimación del modelo Tweedie**
m_tw <- gam(
S_trunc ~ edad_conductor_cat +
edad_vehiculo_cat +
potencia_grupo +
bonus_malus +
combustible +
offset(log(exposicion)),
family = tw(link = "log"),
data = polizas,
method = "REML"
)
summary(m_tw)
##
## Family: Tweedie(p=1.496)
## Link function: log
##
## Formula:
## S_trunc ~ edad_conductor_cat + edad_vehiculo_cat + potencia_grupo +
## bonus_malus + combustible + offset(log(exposicion))
##
## Parametric coefficients:
## Estimate Std. Error t value
## (Intercept) 3.5630992 0.0696660 51.145
## edad_conductor_catAdultos jóvenes -0.2791880 0.0376075 -7.424
## edad_conductor_catAdultos -0.0021799 0.0398750 -0.055
## edad_conductor_catMayores -0.0786176 0.0444586 -1.768
## edad_vehiculo_catMedio 0.0345868 0.0197314 1.753
## edad_vehiculo_catViejo -0.2490611 0.0249864 -9.968
## potencia_grupoBaja -0.2110242 0.0343368 -6.146
## potencia_grupoMedia -0.0920358 0.0329181 -2.796
## potencia_grupoMuy Alta 0.0818471 0.0851122 0.962
## bonus_malus 0.0290061 0.0005877 49.356
## combustibleRegular -0.0847245 0.0182546 -4.641
## Pr(>|t|)
## (Intercept) < 0.0000000000000002 ***
## edad_conductor_catAdultos jóvenes 0.000000000000114 ***
## edad_conductor_catAdultos 0.95640
## edad_conductor_catMayores 0.07701 .
## edad_vehiculo_catMedio 0.07962 .
## edad_vehiculo_catViejo < 0.0000000000000002 ***
## potencia_grupoBaja 0.000000000796443 ***
## potencia_grupoMedia 0.00518 **
## potencia_grupoMuy Alta 0.33623
## bonus_malus < 0.0000000000000002 ***
## combustibleRegular 0.000003463380605 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## R-sq.(adj) = 0.000806 Deviance explained = 3.88%
## -REML = 3.1899e+05 Scale est. = 439.36 n = 677991
Diagnóstico global del ajuste
par(mfrow = c(2,2))
gam.check(m_tw)
##
## Method: REML Optimizer: outer newton
## full convergence after 8 iterations.
## Gradient range [-0.08276832,-0.07467906]
## (score 318991.9 & scale 439.3561).
## Hessian positive definite, eigenvalue range [11029.18,66471.63].
## Model rank = 11 / 11
# Extra: residuos vs ajustados (deviance residuals)
fitted_mu <- predict(m_tw, type = "response")
res_dev <- residuals(m_tw, type = "deviance")
plot(fitted_mu, res_dev,
xlab = "Valores ajustados (μ̂)",
ylab = "Residuos de devianza",
main = "Residuos vs Ajustados")
abline(h = 0, lty = 2)
polizas$mu_freq <- predict(m_nb_final, newdata = polizas, type = "response")
polizas$mu_sev <- predict(mod_sev_final, newdata = polizas, type = "response")
polizas$pp_fs <- polizas$mu_freq * polizas$mu_sev
summary(polizas$pp_fs)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0562 20.8027 57.3075 67.0712 93.4853 6923.8416
quantile(polizas$pp_fs, probs = c(0.5, 0.75, 0.9, 0.95, 0.99))
## 50% 75% 90% 95% 99%
## 57.30750 93.48527 123.45410 165.34326 300.96726
library(ggplot2)
ggplot(polizas, aes(x = log(pp_fs))) +
geom_histogram(bins = 50, fill = "steelblue", color = "white") +
labs(
title = "Distribución de la prima pura estimada (Frecuencia–Severidad)",
x = "log(Prima pura estimada)",
y = "Frecuencia"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 12),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
df_conc <- polizas %>%
filter(is.finite(pp_fs)) %>%
arrange(desc(pp_fs)) %>%
mutate(
prop_polizas = row_number() / n(),
prop_coste = cumsum(pp_fs) / sum(pp_fs)
)
ggplot(df_conc, aes(x = prop_polizas, y = prop_coste)) +
geom_line(linewidth = 1, color = "black") +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "steelblue") +
labs(
title = "Concentración del coste esperado",
x = "Proporción acumulada de pólizas (ordenadas de mayor a menor)",
y = "Proporción acumulada del coste esperado"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 12),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
polizas <- droplevels(polizas)
mod_tw_final <- gam(
S_trunc ~ edad_conductor_cat +
edad_vehiculo_cat +
potencia_grupo +
bonus_malus +
combustible +
offset(log(exposicion)),
family = tw(link = "log"),
data = polizas,
method = "REML"
)
polizas$pp_tw <- predict(mod_tw_final, newdata = polizas, type = "response")
summary(polizas$pp_tw)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1803 25.6107 70.8986 82.1729 114.9888 6148.0818
ggplot(polizas, aes(x = log(pp_tw))) +
geom_histogram(bins = 50, fill = "steelblue", color = "white") +
labs(
title = "Distribución de la prima pura estimada (Tweedie)",
x = "log(Prima pura estimada)",
y = "Frecuencia"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 12),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8)
)
mean(polizas$pp_fs)
## [1] 67.07117
mean(polizas$pp_tw)
## [1] 82.1729
sum(polizas$pp_fs)
## [1] 45473647
sum(polizas$pp_tw)
## [1] 55712485
tabla_cuantiles <- rbind(
FS = quantile(polizas$pp_fs, probs = c(0.5, 0.75, 0.9, 0.95, 0.99)),
TW = quantile(polizas$pp_tw, probs = c(0.5, 0.75, 0.9, 0.95, 0.99))
)
tabla_cuantiles
## 50% 75% 90% 95% 99%
## FS 57.30750 93.48527 123.4541 165.3433 300.9673
## TW 70.89863 114.98877 146.2145 203.3443 373.2328
cor(polizas$pp_fs, polizas$pp_tw)
## [1] 0.9874853
polizas$pp_emp <- polizas$importe_medio_siniestro / polizas$exposicion
comp_agregada <- polizas %>%
summarise(
n = n(),
O_total = sum(pp_emp, na.rm = TRUE),
E_FS_total = sum(pp_fs, na.rm = TRUE),
E_TW_total = sum(pp_tw, na.rm = TRUE),
O_media = mean(pp_emp, na.rm = TRUE),
E_FS_media = mean(pp_fs, na.rm = TRUE),
E_TW_media = mean(pp_tw, na.rm = TRUE),
O_mediana = median(pp_emp, na.rm = TRUE),
E_FS_mediana = median(pp_fs, na.rm = TRUE),
E_TW_mediana = median(pp_tw, na.rm = TRUE)
)
comp_agregada
## n O_total E_FS_total E_TW_total O_media E_FS_media E_TW_media
## 1 677991 243491667 45473647 55712485 359.137 67.07117 82.1729
## O_mediana E_FS_mediana E_TW_mediana
## 1 0 57.3075 70.89863
cal_fs <- polizas %>%
filter(is.finite(pp_fs), is.finite(pp_emp)) %>%
mutate(decil = ntile(pp_fs, 10)) %>%
group_by(decil) %>%
summarise(
n = n(),
O = mean(pp_emp, na.rm = TRUE),
E = mean(pp_fs, na.rm = TRUE),
OE = O / E
)
cal_tw <- polizas %>%
filter(is.finite(pp_tw), is.finite(pp_emp)) %>%
mutate(decil = ntile(pp_tw, 10)) %>%
group_by(decil) %>%
summarise(
n = n(),
O = mean(pp_emp, na.rm = TRUE),
E = mean(pp_tw, na.rm = TRUE),
OE = O / E
)
ggplot(cal_fs, aes(x = decil, y = OE)) +
geom_line() + geom_point() +
geom_hline(yintercept = 1, linetype = "dashed") +
labs(title = "Calibración por deciles (FS)", x = "Decil (PP_FS)", y = "O/E") +
theme_minimal()
ggplot(cal_tw, aes(x = decil, y = OE)) +
geom_line() + geom_point() +
geom_hline(yintercept = 1, linetype = "dashed") +
labs(title = "Calibración por deciles (Tweedie)", x = "Decil (PP_TW)", y = "O/E") +
theme_minimal()
write.xlsx(polizas, "base_datos.xlsx")